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Turbulent transport in the solar nebula 

By K. W. THOMPSON 

1. Problem Description 

It is likely that turbulence played a major role in the evolution of the so- 
lar nebula, which is the flattened disk of dust and gas out of which our solar 
system formed. Relevant turbulent processes include the transport of angular 
momentum, mass, and heat, which were critically important to the formation 
of the solar system. This research will break new ground in the modelling of 
compressible turbulence and its effects on the evolution of the solar nebula. The 
computational techniques which have been developed should be of interest to re- 
searchers studying other astrophysical disk systems (e.g. active galactic nuclei), 
as well as turbulence modelers outside the astrophysics community. 

2. Objectives/Milestones 

The work currently being performed system is closely connected with that of 
Cabot, Hubickyj, and Pollack (hereafter CHP), who have been using the turbu- 
lent channel flow code of Kim, Moin, & Moser (1987) to investigate incompress- 
ible turbulence in the solar nebula. The turbulence simulation code developed 
for this project will ultimately replace the central elements of the channel flow 
code, while incorporating many of the useful analysis tools developed for the 
latter. 

The objectives of this project fall into two categories: 

1. To advance the understanding of the role of compressible turbulence in the 
evolution of the solar nebula by carrying out a series of numerical experiments 
to evaluate the Reynolds stress tensor, turbulent heat transfer rate, turbulent 
dissipation rate, and turbulent energy spectrum for conditions relevant to the 
solar nebula; 

2. To advance the state of the art in turbulence simulation techniques by 
exploring new computational methods which are expected to be both efficient 
and flexible. 

The field of solar nebula physics is one of much current interest, as seen, for 
example, in the contributions to Protostars and Planets II (Black & Matthews 
1985). The field of turbulence is also one of much current interest. 

2.1 Scientific Objectives 

The first objective is scientific and largely astrophysical in nature. Turbulent 
processes are believed to have been important in the evolution of the solar nebula, 
which was a rarefied disk of gas and dust, out of which the planets, asteroids, and 
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comets formed. This disk circled the Sun during and shortly after its formation. 
The nebula is in turn thought to have formed out of the contraction of a much 
larger and even more diffuse molecular cloud. The nebula’s central star, our 
Sun, also formed out of the molecular cloud matter, and it is the Sun’s gravity 
which held the nebula together and kept it from flying apart. The combination 
of the Sun’s gravitational field, the initial angular momentum of the molecular 
cloud (retained by the nebular material), and radiative cooling is believed to have 
confined the solar nebula to a thin disk, rather than a cloud. (The formation and 
evolution of protostellar disks continue to occur in active star formation regions 
in our galaxy, and this research is relevant to nebular problems in general, not 
only to our solar system.) 

Anticipated sources of turbulence include thermal convection (driven by tem- 
perature gradients), vertical shear in the angular velocity (due to the nonuniform 
deviations from the central plane orbital motion of the gas with altitude), and 
the large velocity differential between the rotating disk matter and the infalling 
(non-Keplerian) molecular cloud material. Turbulent convection is expected to 
be strongly dependent on the compressibility of the flow; hence, the study of 
compressible turbulent flow is crucial to this project. 

The effects of turbulence on the solar nebula are many. Turbulent dissipation 
causes a loss of orbital energy and, therefore, an inward flow of mass toward the 
central star. The inflow of mass must be accompanied by an outflow of specific 
angular momentum in order to conserve total angular momentum. The inward 
transport of mass feeds most of the nebular matter to the Sun, while the outward 
transport of angular momentum causes the nebula to expand. Turbulence will 
also transport heat, both radially (parallel to the central disk plane) and verti- 
cally (in the perpendicular direction). Turbulent dissipation, mentioned above, 
is expected to be a major (perhaps dominant) source of heat for the nebula. 
Turbulent dissipation converts kinetic energy to heat, which is then radiated 
away to space, reducing the total energy of the system and making the nebula 
more tightly bound to the Sun. 

The greatest obstacle to the accurate modeling of the solar nebula is the fact 
that the length scales on which viscous dissipation takes place (and on which 
the turbulent kinetic energy is turned into heat) is many orders of magnitude 
smaller than the size of the disk. Consequently, it is not currently possible to 
create a single computational model which accurately simulates both the large 
scale structure and the small scale turbulent dissipative processes. The objective 
of this project is to simulate numerically the turbulent processes on a small scale 
and obtain a parameterization of these processes which may be used in other 
attempts to model the large scale evolution of the nebula. Much of this work 
is being done by other researchers (CHP), but their efforts must be considered 
in the developments described below, as the results of their research will be 
incorporated into this model. 

Several issues need to be addressed in order to produce a realistic model for 



Turbulent transport in the solar nebula 


109 


the solar nebula. The most obvious of these is the Sun’s gravity, which forces 
the gas to follow nearly Keplerian velocity profiles in its rotation about the Sun. 
Radial shear is, therefore, produced in the flow, as the orbital angular velocity 
decreases with distance from the Sun. 1 

The Keplerian flow velocities are highly supersonic in the rest frame of the 
central star. A direct simulation of the flow in the star’s rest frame is impractical, 
as the flow velocities dictate unworkably small time steps. A better approach 
is to work in a coordinate system which is comoving with the average flow in 
the model volume, so that the velocities are subsonic in the model’s coordinate 
system, and the time steps are more reasonable. 

A coordinate transformation based on the work of Rogallo (1981) will be used 
to represent the radial shear in a form which permits the use of periodic boundary 
conditions in the radial direction. The Rogallo transformation eliminates the 
need to devise boundary conditions which properly advect turbulent flow in and 
out of the radial boundaries. 

A second type of shear exists due to the matter falling in from the collapsing 
molecular cloud onto the disk. The radial and angular velocity components 
of the infalling material will be quite different from those of the rotating disk. 
While the vertical velocity component of the infalling matter will be reduced by 
the accretion shock at the cloud-disk interface, it will remain nonzero, providing 
a downward directed ram pressure on the disk material. In addition, the radial 
and singular velocity components of the infalling matter will not be affected by 
the shock and will create a large velocity shear relative to the disk material. 
This is an important question, and new techniques may have to be developed to 
study it. The effects of these factors on nebular turbulence could be important. 

A subtler gravitational effect is the vertical variation of gravity within the 
disk. The disk is very thin compared to its radial dimensions; therefore, a 
fluid element a distance z above the central disk plane experiences a downward- 
directed gravitational force (due to the central star) which is proportional to 
z. This linear variation of gravity is expected to have a significant effect on 
the convective flow. Convection in the Earth’s atmosphere takes place in an 
altitude range over which the gravitational acceleration is essentially constant. 
The nebular problem has a variable acceleration with altitude, which may give 
rise to flow patterns not seen in constant gravity environments and is expected 
to affect turbulent transport. 

Another major issue is the compressibility of the fluid. The computational 
technique of Kim et al, as currently implemented by CHP assumes an incom- 
pressible fluid, but the solar nebula is strongly stratified and is, therefore, ex- 
pected to behave like a compressible gas. Thus it is necessary to use the com- 
pressible fluid equations. Most of the work on turbulent flow in the past has 

1 Note: the radial shear by itself is stable and will not directly cause turbulence. Turbulence 
must be induced by other effects such as viscous heating and vertical shear. 
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concentrated on incompressible flows. 

The channel flow code assumes flow between two parallel plates, for which the 
fluid velocities go to zero at the plates. The no-slip boundary condition makes 
sense for flow in a channel but is not valid for flow in the solar nebula. It is 
desirable to replace these no-slip wall boundary conditions with conditions which 
are representative of the solar nebula. The nonreflecting boundary conditions 
of Thompson (1987a) are being used for this purpose. These conditions allow 
material to flow across the boundaries and allow outward propagating waves to 
leave the computational domain without generating reflections. 

2.2 Improvements to Turbulence Simulation Techniques 

The second category of objectives, which may have applications outside of 
astrophysics (e.g. the aeronautical community), focuses on the improvement of 
turbulence simulation techniques in general. 

The computer code currently in use is very flexible and is able to adjust 
to a wide variety of boundary conditions and computational techniques. It is 
designed to be easily modified so that new features may be added and differ- 
ent computational techniques incorporated as needed. For example, the spa- 
tial derivative methods can be modified at will (e.g. from finite difference to 
pseudospectral) by altering only that section of the code which computes these 
derivatives. 

The techniques in use are described in section 3. 

2.3 Milestones Passed 

The work performed so far has been concerned with analyzing alternative 
computational methods, developing a suitable computing strategy, and validat- 
ing the computer simulation. Early design and testing focused on simple linear 
and nonlinear wave simulations, for which error analysis and convergence testing 
could be performed easily. The later stages have involved the full set of compress- 
ible fluid equations. Successful tests have validated the design of the computer 
code and are described below. One of these tests has had the added benefit of 
generating a significant new research project in its own right, as described below. 

The major program tests are: 

1. Rarefaction wave. 

2. Oscillating atmospheric column. 

3. Hydrodynamic escape. 

4. Unstable shear flow (Kelvin-Helmholtz instability). 

2.3.1 Rarefaction Wave 

The rarefaction wave is the well known solution to the adiabatic expansion 
of a uniform gas into a vacuum in the absence of gravity. The exact solution 
to this problem has been given in many places, (e.g. Landau & Lifshitz 1959, 
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Thompson 1986). Time dependent numerical simulations of the continuous part 
of the rarefaction wave were highly accurate and exhibited the proper fourth 
order convergence with grid refinement. 

2.3.2 Oscillating Atmospheric Column 

This problem describes a column of atmosphere in a constant gravitational 
field. The density and pressure drop off exponentially with altitude z , as e~ z ! h , 
where h is the vertical scale height. This is a one dimensional problem, as the 
solution does not depend on x or y. The initial state is isothermal and is stable. 

Just as in an organ pipe, one can set up standing waves at certain discrete fre- 
quencies. I have set up an adiabatic perturbation to excite the lowest frequency 
mode of the column for problems ranging from 1 to 50 vertical scale heights 
in extent. In all cases the time dependent solution accurately reproduced the 
expected oscillatory mode and displayed the correct fourth order convergence. 
The exponential dependence of the unperturbed state makes this problem an 
obvious choice for solution by the methods of section 3.1. 

2.3.3 Hydrodynamic Escape 

This problem is also one dimensional but spherically symmetric. It describes 
the escape of a planetary atmosphere into space. Gravity obeys the inverse 
square law, and the flow velocity goes smoothly from highly subsonic near the 
planet’s surface to highly supersonic at large distances. Steady state solutions 
have been known analytically for some time (Bondi 1952). It is an interesting 
fact that the same equations describe the solar wind and the accretion of gas 
onto a protostar embedded in a spherically symmetric cloud. 

Although the analytic solution is known, this problem presents special prob- 
lems for numerical techniques (Zahnle 1988) . I have found it to be a challenging 
problem (and a useful one for uncovering some subtle boundary condition prob- 
lems) but can now produce stable and accurate solutions. The problem contains 
enormous pressure and density gradients near the lower boundary, and the den- 
sity and pressure vary by several orders of magnitude in the model volume. This 
problem is rendered tractable by the logarithmic techniques described in section 
3.1. Previous attempts to solve this problem numerically with more standard 
methods have failed. 

The more general problem of hydrodynamic escape of planetary atmospheres 
which are affected by thermal conductivity and solar heating is of considerable 
interest to researchers in the field. The successful solution of the more basic 
problem will serve as a starting point for further research on this subject by 
Zahnle, my co-investigators, and myself. 

2.3.4 Unstable Shear Flow 

This is a well known problem in the atmospheric sciences, and one which is 
two dimensional in nature. We have an initially isothermal atmosphere whose 
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density and pressure vary exponentially with altitude z as e~ x ! h (rectangular 
coordinates). The atmosphere contains a shear flow with velocity components 
(U(z), 0,0), where U(z) = U 0 tanh(z/d), d < h, and Uq is subsonic. We perturb 
the flow field by adding small velocity perturbations u(x, z) and w(x,z) so that 
the perturbed velocity is ( U + u,0, u;). Perturbation theory (Chandrasekhar 
1961) shows that the stability of the flow is determined by the Richardson number 
Ri = gcP/hU^. A Kelvin-Helmholtz instability arises if Ri < 1/4, while the flow 
is stable if Ri >1/4. 

Calculations performed with Richardson numbers of 0.05 and 0.15 do indeed 
give rise to unstable flow. These calculations show that the initially small pertur- 
bations grow and create a vortex. Similar calculations with Richardson numbers 
of 0.35 and 1.25 are stable and do not lead to vortex formation. These results 
are consistent with the predictions of linear stability analysis. 

3. Technical Approach 

The equations to be solved are the compressible Navier Stokes equations. 
They describe the time evolution of a compressible fluid in three dimensions and 
incorporate the effects of a variable gravitational field, radial shear, molecular 
viscosity, thermal conductivity, and internal heat sources. The solution process 
is made challenging by the large variation in density and pressure which occurs in 
this problem. A new technique has been developed to handle the range problem 
as described below. 

The purpose of this calculation is to study relatively small scale turbulence 
in order to characterize the effects of turbulence on the overall flow without 
resorting to ad hoc turbulent viscosity approximations. Therefore, the actual 
volume to be simulated consists of a small “box” of nebular material in the disk 
of relatively small radial and angular extent as seen from the Sun. 2 The box 
includes one to several scale heights of pressure variation in the direction per- 
pendicular to the disk midplane. Periodic boundary conditions are used at the 
boundaries facing “sideways” into the disk material (i.e. the (f> boundaries) . The 
r boundaries make use of Rogallo’s method (Rogallo 1981) to specify periodic 
boundary conditions. The d (vertical) boundaries require non-trivial boundary 
conditions and are currently handled by the nonreflecting conditions described 
by Thompson (1987a). 

3.1 The Range Problem 

The density and pressure in the nebula fall off very rapidly with altitude above 
the nebula midplane, roughly as , where h is a scale height. 3 Thus the 

density and pressure may vary by orders of magnitude throughout the model 
volume, which poses a difficult challenge to conventional numerical techniques. 

The spherical coordinate system is assumed for this discussion. 

o 

The unusually rapid falloff is due to the variable gravity, which increases with altitude. 
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I have developed a technique to deal with scalar fields which vary over such a 
large range by using the logarithms of density and pressure instead of the basic 
quantities themselves. For example, the one dimensional continuity equation, 


dp , d 


is replaced by 


dlnp 3 In p du 

HT + u ~aT + ai = 0. 


and similarly for the energy equation. Thus an exponential dependence of p 
on x is reduced to a polynomial dependence of In p on x, which the numerical 
scheme can cope with more readily. All derivatives of density, pressure, and 
thermal energy may be replaced by their logarithmic forms as above, leading 
to a more tractable system. An added benefit is that the density and pressure 
cannot become negative with this approach. 


3.1.1 The Numerical Techniques 

The current code uses standard fourth order finite difference formulas to eval- 
uate spatial derivatives in all directions (see, for example, Thompson 1987b). 4 

Time integration is handled by the usual fourth order explicit Runge-Kutta 
method, which is simple to implement and has excellent stability properties. 
This combination of spatial derivative and time integration methods provides 
global fourth order convergence. The fourth order convergence rate has been 
verified repeatedly on a large number of test runs and implies that this code 
requires significantly fewer grid points than a second order code to achieve a 
given level of accuracy. 

The choice of an explicit method over an implicit method stems from the need 
to resolve the smallest features present in the flow. At the smallest length scales, 
viscosity dominates the evolution of the flow. Since we need to simulate the 
dissipation of kinetic energy to heat accurately at these scales, the grid spacing 
and time steps necessary are set by the properties of the flow and are the same 
whether explicit or implicit methods are used. The optimal grid spacing and 
time step are those for which the propagation and viscous Courant numbers 
are equal. Consequently, the simpler (and faster) explicit approach has been 
selected. 


4 Compact difference formulas of the compact (Pade) type (Rubin & Koshla 1977) have been 
evaluated but have not been found to offer significant improvements over the standard fourth 
order formulas. 
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